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Abstract 

The putative recent indication of an unidentified 3.55 keV X-ray line in certain astrophysical 
sources is taken as a motivation for an improved theoretical computation of the cosmological 
abundance of 7.1 keV sterile neutrinos. If the line is interpreted as resulting from the decay 
of Warm Dark Matter, the mass and mixing angle of the sterile neutrino are known. Our 
computation then permits for a determination of the lepton asymmetry that is needed for pro¬ 
ducing the correct abundance via the Shi-Fuller mechanism, as well as for an estimate of the 
non-equilibrium spectrum of the sterile neutrinos. The latter plays a role in structure forma¬ 
tion simulations. Results are presented for different flavour structures of the neutrino Yukawa 
couplings and for different types of pre-existing lepton asymmetries, accounting properly for 
the charge neutrality of the plasma and incorporating approximately hadronic contributions. 
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1. Introduction 


Despite considerable phenomenological success, the Standard Model of particle physics suffers 
from a number of shortcomings: it can explain neither the presence of neutrino oscillations, 
nor of dark matter, nor of a cosmological baryon asymmetry. It is remarkable that all of 
these problems can in principle be cured by a simple enlargement of the theory through three 
generations of right-handed neutrinos [1]”[4], without changing the underlying theoretical 
principles of gauge invariance and renormalizability. Despite its simplicity it remains unclear, 
of course, whether nature makes use of this possibility. 

In the present paper we are concerned with the dark matter aspect (for reviews see, e.g., 
refs. [5, 6]). For the parameter values that are phenomenologically viable, the dark matter 
sterile neutrinos do not contribute to the two observed active neutrino mass differences, so 
that the dark matter aspect is partly decoupled from the neutrino oscillation and baryon 
asymmetry aspects. The decoupling is not complete, however: it turns out that the sterile 
neutrino dark matter scenario is tightly constrained [7], and only works if the dynamics 
responsible for generating a baryon asymmetry also generates a lepton asymmetry much 
larger than the baryon asymmetry [8], which subsequently boosts sterile neutrino dark matter 
production through an efficient resonant mechanism first proposed by Shi and Fuller [9]. 
Despite the tight constraints, indications of a possible observation [10, 11] demand us to 
take this scenario seriously. For our purposes, the complicated dynamics of ref. [8] simply 
amounts to the fact that the initial state at a temperature of a few GeV possesses certain 
lepton asymmetries. 

The goal of the present paper is to rehne the analysis of ref. [7] in a number of ways, both 
theoretically and as far as the numerical solution of the rate equations is concerned. The 
resulting spectra could be used as starting points in structure formation simulations (cf. e.g. 
refs. [12]-[14] and references therein). 

Our plan is the following. In sec. 2 we review and rehne the theoretical description of 
sterile neutrino dark matter production from a thermalized Standard Model plasma with pre¬ 
existing lepton asymmetries. The practical implementation of the corresponding equations 
and a strategy for their solution are discussed in sec. 3. Numerical solutions are presented in 
sec. 4, and we conclude with a discussion in sec. 5. 


2. General derivation of the rate equations 

Our hrst goal is to derive a closed set of rate equations for the sterile neutrino distribution 
function and for lepton number densities,^ valid both near and far from equilibrium. In ref. [7] 
such equations were postulated but the argument involved phenomenological input, in order 

^By “number densities” we mean asymmetries, i.e. particles minus antiparticles. 
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to account for two types of “back reaction”. Our derivation yields an outcome that differs 
slightly from ref. [7]. Moreover, the influence of electric charge neutrality of the plasma on 
the relation between lepton number densities and lepton chemical potentials was not properly 
accounted for in ref. [7]. Finally, in order to obtain a maximal effect, the lepton asymmetries 
of the different flavours were treated as equilibrated in ref. [7], even though flavour equilibrium 
(through active neutrino oscillations) is expected to be reached only at temperatures below 
10 MeV or so [15, 16]. We eliminate all of these simplifications in the present paper. 

We start by defining, in sec. 2.1, the quantities appearing in the equations. In sec. 2.2 an 
evolution equation is obtained for the sterile neutrino distribution function. In sec. 2.3 the 
same is achieved for lepton number densities, with a right-hand side expressed in terms of 
lepton chemical potentials and sterile neutrino distribution functions. The system is com¬ 
pleted in sec. 2.4, where we relate lepton chemical potentials to lepton densities, taking into 
account electric charge neutrality of the Standard Model plasma. 


2.1. Definitions 


We consider a system consisting of right-handed (sterile) neutrinos and Standard Model (SM) 
particles. The SM particles are assumed to be in thermal equilibrium at a temperature T. 
The initial state is characterized by non-zero lepton chemical potentials, associated with 
different lepton flavours. (At low temperatures T < 10 MeV the lepton flavours equilibrate 
through active neutrino oscillations, so that there is only a single lepton chemical potential, 
denoted by /U^.) The initial state may also contain an ensemble of sterile neutrinos. The 
two sectors communicate through Yukawa interactions, parametrized by neutrino Yukawa 
couplings. As a result, the distribution function of the sterile neutrinos evolves towards its 
equilibrium form, and the lepton densities decrease, because lepton number is violated in 
the presence of neutrino Yukawa interactions and Majorana masses. (However, in practice 
neither process gets completed within the lifetime of the Universe.) 

Let us denote by h the matrix of neutrino Yukawa couplings. We work out a set of rate 
equations to 0{h?). This means that, as soon as we have factorized a coefficient of 0{h?), 
the sterile neutrinos can be treated as mass eigenstates and free particles, with masses given 
by Majorana masses. Sterile neutrinos can then be represented by on-shell field operators. 


Nj{X) 

N,{X) 


d^k 


d^k 


+ “/ko-^/kcr® 




/ , “/ko-^/ko-e + “/ko-^/ko-e 


( 2 . 1 ) 

( 2 . 2 ) 


where Ej = and JCj - X = Ejt — • X. The index I enumerates the sterile 

neutrino “flavours”, and M/ is their mass, which we assume to be real and positive. The 
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creation and annihilation operators satisfy the anticommutation relations ~ 

— p), consistent with {iVj(t, x), ivj(t, y)} = — y). The on-shell spinors 

u,v are normalized in a usual way, for instance ^^/ko-^/kr ~ ~ 

'^a^ika^ika = ^ + Mj and ^/kcr^/kcr — f' ~ Majorana nature of the spinors 

requires that u = Cv^, v = Cu^, where C is the charge conjugation matrix. 

Now, let us define the ensemble occupied by the sterile neutrinos. Their distribution is 
described by a density matrix, denoted by pj^. We take a rather general ansatz for pj^, 
assumed however to be diagonal in flavour, momentum, and spin indices: 

PN = Z-^expn M/k«ka«/ka) > (2-3) 

I,a 

where the function is left unspecified except for being spin-independent, Zj^ takes care 
of overall normalization, and fj^ = f d^k/(27r)^. The density matrix has also been assumed 
to be x-independent. Note that even though the function pj^. bears some resemblance to a 
chemical potential, it is not identical to one (Majorana fermions are their own antiparticles 
and no chemical potential can be assigned to them). 

The property originating from pj^ that we need in practice is a phase space distribution 
function, denoted by fj^^, which can be defined as 

Tl' i^ma^JprPN) = (k “ P) //k • (2-4) 


Setting the indices equal, this amounts to 

fji, = (4ka«JkaPjv) ^ (2.5) 

where V denotes the volume. We remark that the normalization in eq. (2.5) differs from that 
in ref. [7] by a factor (27r)^, and is such that the total number density of sterile neutrinos, 
summed over the two spin states, reads 


nj = 



( 2 . 6 ) 


Let us then turn to the Standard Model (SM) part. For its density matrix we adopt the 
ansatz 


PsM — -^SM 


rp 'y ^ Pa^a 


(2.7) 


where B is the baryon number operator (at T <C 160 GeV baryon and lepton numbers are 
separately conserved within the SM). The lepton number operator associated with flavour a 
reads 


= 


ea7o(aL + 


( 2 . 8 ) 
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Here Ca denotes a charged lepton of flavour a, and 0 ^ = (1 — 75 ) 72 , a^= (1 + 75)72 are chiral 
projectors. The left-handed doublet is denoted by Ga)'^■ 

The interaction between SM degrees of freedom and sterile neutrinos is described by the 
interaction Hamiltonian 


— 


I ,a 


/ {^lhlaja+ jah}a^l) 


(2.9) 


where /i is a 3 x 3 Yukawa matrix with complex elements /ij^. The gauge-invariant operators 
to which the sterile neutrinos couple are 

3a = ‘P^aJa > 3a = > ( 2 - 10 ) 


where cj) = i(T 2 (p* is the conjugate Higgs doublet. For the moment the only property that is 
needed from SM dynamics is the spectral function corresponding to these operators, 

PabilC) = {j,(Y), i5(0)}) . (2.11) 

In practice we assume this function to be flavour-diagonal, i.e. oc however the weight is 
flavour-dependent because charged lepton masses play an important role.^ 

For describing the dynamics of the coupled system, we start from an “instantaneous” initial 
state, 

p{^) = Psm®Pn, ( 2 - 12 ) 

where the SM part is from eq. (2.7) and the sterile neutrino part is from eq. (2.3). In an 
interaction picture, the density matrix evolves as 

p,{t)=p{Q)-i fdt' [HAt'),m] - fdt' fdt" +0{h^) , (2.13) 

3o 3o 3o 

where is the Hamiltonian corresponding to eq. (2.9) in the interaction picture. This evo¬ 
lution causes both the SM and sterile neutrino density matrices to evolve. For the Standard 
Model this change is almost negligible whereas for sterile neutrinos it is of 0(1). 

In order to obtain differential equations for the set of observables to be defined, the time t 
is considered large compared with time scales of the SM plasma, t S> l/{a‘^T), where a is a 
generic fine-structure constant. This guarantees that decoherence takes place in the SM part 
of the Hilbert space. At the same time t should be small enough to avoid secular terms. In 
practice, through an appropriate choice of and p^, we can arrange things in a way that 
the limit t —)■ oo can be taken and equilibration is correctly accounted for without secular 
terms (see below). An implicit assumption we make is that since the sterile neutrinos are 

^Note that in the range of temperatures relevant for us, T 7> 10 MeV, active neutrino flavour oscillations 
are not fast enough to play a role [15, 16], and we can work directly in the interaction basis. 
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being produced from a statistical plasma, they are decoherent, and their density matrix is 
assumed to retain the diagonal form in eq. (2.3); this assumption would surely be violated 
by terms of higher order in h. However, given that we do find a system evolving towards 
equilibrium, the ansatz of eq. (2.3) can be considered self-consistent. 

2.2. Sterile neutrino production rate 

As a first ingredient, we determine the production rate of sterile neutrinos from an asymmetric 
plasma. The computation can be carried out with the formalism of ref. [17].^ From eq. (2.13) 
the density matrix of the system evolves as 

/5j(t) = (odd powers of /i) — f dt' [.H;(t^),/5(0)]] +0{h‘^) , (2-14) 

Jo 

where terms containing odd powers of h have been suppressed because they are projected 
out later on. The density matrix of the initial state is given by eq. (2.12). 

The observable we are interested in corresponds to the time derivative of the expectation 
value of eq. (2.5), normalized like in eq. (2.4) (and summed over spin states): 

2//k = ^ (4ka«/kaA) • (2-15) 

CF 

Inserting eq. (2.14), we are faced with 4-point functions of the creation and annihilation 
operators, evaluated in an ensemble defined by pj^. Given that the density matrix is assumed 
diagonal in flavour, momentum and spin indices and has an effectively Gaussian appearance. 


the 4-point functions can be reduced to 2-point functions;^ 

Tt i^rq^kp ^rp^kq)frfk ’ (2-16) 

Tr (dJdpa^cigPjy) Sj,p6j^gf^fp. /^) , (2.17) 

Tr (cLl,Clp(lqCLy.p jq) i^rp^kq ^rq^kp)fri^ fk) ’ (2-18) 

Tr {cipCl\,ClqCly.Pfq^ (5j,p(5^g(l /fc T frfk) T ^rq^kpfri^ /fc) > (2-19) 

Tr {apdlaldgP^) = (5^p4q(l - fr)fk - ^rq^pfrC^ - fk) > (2-20) 


where the indices incorporate all dependences. Subsequently a somewhat tedious analysis, 
tracing the steps outlined in ref. [17], shows that all terms quadratic in the distribution 
functions drop out. The final result reads 

2//k(^) = “ ^a) “ flM] Tr [If Paa{Kf) ^r] 

a ^ 

+ [nY{Ei + Pa) - //k(*)] Tr [If Paai-E) Or] } > (2-21) 

®It could also be carried out with the (non-equilibrium) Schwinger-Keldysh formalism. In order to ob¬ 
tain kinetic equations, one may Wigner-transform Schwinger-Dyson equations for the forward and backward 
Wightman propagators of the sterile neutrinos (see e.g. ref. [18] for the generic formalism). 

^Depending on how the computation is organized, some of these relations may not be needed. 
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where np is the Fermi distribution. On the right-hand side we have replaced //^(O) with 
with p^{t) in eq. (2.14). This is consistent with our goal of an 0{h?) 
determination of fj^, since the difference between /Oi(0) and Pi{t) is of higher order in h. 
As was anticipated at the end of the previous section, this prescription removes any secular 
terms from the evolution. Finally, we remark that, in the limit fj^ 0, the result agrees 
with ref. [7], whereas in equilibrium, corresponding to = 0, fj^ = np{Ej), the production 
rate vanishes as must be the case. 


2.3. Lepton number washout rate 

The lepton number defined by eq. (2.8) is not conserved because of the interactions in eq. (2.9). 
The operator equation of motion reads [19] 


ia = lj(^NhT^j -jn^N) , 


( 2 . 22 ) 


where {T°')ij = ^ai^aj- propose to evaluate the expectation value of this operator (taken 
in the interaction picture) in the time-dependent ensemble described by eq. (2.13). This time 
the leading contribution comes from the term linear in = Lg/V): 


Tr 


n„ = 


lim — 

L 

V 

t^oo 

lim — / 

t—>-oo Y 

■^x,y 


LaPlit) 


dt' 


{NhT^j -jT^h^N){X), {Nhj + jh^N){y)\ ) + 0{h^) , (2.23) 


where X = (t,x), y = (t',y), and (...) = Tr [(...)/5(0)]. The SM part of this expectation 
value can be expressed in terms of Wightman correlators, which in turn can be expressed in 
terms of the spectral function in eq. (2.11). For the sterile neutrino part, we can insert the 
field operators from eqs. (2.1), (2.2), and eliminate the annihilation and creation operators 
through eq. (2.4). A tedious but straightforward analysis yields 


= ^^^^^{[//k(^) -™f(-Ej -/r„)]Tr [^p^„(/C)aJ 

+ [n^Ej + pg) - TV pgg{-lC) a J } + 0{h^) , (2.24) 

where we have employed the prescription described below eq. (2.21). In the limit —)■ 0, 
the result corresponds to eq. (2.32) of ref. [7]. On the other hand in equilibrium, i.e. Pg = 0, 
= np{Ej), the washout rate vanishes as must be the case. 

An interesting crosscheck of eq. (2.24) can be obtained by putting sterile neutrinos in 
equilibrium {fj^^ —)■ np{Ej)) and by expanding lepton number densities to first order around 
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equilibrium. In the limit of small we can re-express the through the lepton densities 
through a susceptibility matrix 1^6=0’ Then we get 

< = -lah^b + n^) , (2.25) 


describing how lepton asymmetries disappear (or are “washed out”) near equilibrium. The 
coefficients read 


Tafe 


\hla 


<{E,] 

Er 


Tr 


[jf [pM + Paai-E)] a^} , (2.26) 


which after an appropriate adjustment of conventions agrees with eq. (4.7) / (29) of ref. [19]. 


2.4. Relation of lepton densities and chemical potentials 


In order to close the set of equations, we need to determine the relations between lepton 
densities and chemical potentials. Even though the densities are separately conserved 
within the Standard Model with h = 0, their fluctuations are correlated. The reason is 
that, because of charge neutrality, an excess of electrons over positrons is compensated for 
by an excess of antimuons over muons. This implies that the relation of pa and Ua is non¬ 
diagonal. In the present section we work it out to leading order in Standard Model couplings, 
at T < 1 GeV. (Recently, similar computations have been extended up to higher orders in 
Standard Model couplings at T> 160 GeV [19, 20].) 

The desired relations can most conveniently be obtained by first computing the pressure 
(i.e. minus the grand canonical free energy density, p = —Vt/V) as a function of chemical 
potentials associated with all conserved charges. Apart from baryon and lepton numbers, 
chemical potentials need to be assigned to gauge charges, in our case the electromagnetic 
U(l)gm charge (= pq) and the weak SU(2 )l charge {pz) [21]- The chemical potentials assigned 
to gauge charges correspond to the fact that the zero components of the associated gauge 
helds can develop an expectation value. In our case, the weak gauge bosons can be omitted, 
because their effective potential has a large tree-level term ~ p\v‘^, which implies that Pz 
is suppressed by with respect to pq. Therefore, only pa,PQ, and pq need to be 

included. Given that the chemical potentials are very small compared with the temperature 
{pa < 0.02T), it is sufficient to determine p up to 0{p^) in the chemical potentials. 

Omitting exponentially suppressed contributions from the gauge bosons and from the 
top quark, the pressure can be expressed as 


P{Pa,PB,PQ) 

^P{Pa,PB,PQ) 


= p{0) + Ap{p^,PB,Pq) + O{p‘^), 


(2.27) 


= Nr 


Pb Pq 


Xirrii 

i=d,s,b 


+ X/ (Pa-PQ)‘^X{'>^ea) +0{a^) 


(2.28) 


a=e,fi,T 
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where A^c = 3 is a book-keeping variable for hadronic effects, and x is a (diagonal) “suscep¬ 
tibility”. In the free limit the susceptibility reads 


X(m) = 2 / [- 


[ i-nUE)] 

Jp 


d^p 2nF(il^)[l — np{E)] 
(27r)3 f 


, E=y/^ 




(2.29) 


For vanishing mass this evaluates to x(0) = whereas for a non-vanishing mass it can 

be expressed as 


2 oo 




(2.30) 


n=l 


where K 2 is a modified Bessel function. (Radiative corrections to diagonal quark susceptibili¬ 
ties have been computed up to a high order in the massless limit [22], and the same quantities 
can also be measured on the lattice, cf. e.g. refs. [23, 24] and references therein.) 

Given eq. (2.28), is eliminated in favour of the baryon density through = dp/dfi^, 
yielding 


2Nr 


Ub = 


9 1 


Xucil^B + 2liQ) + Xdsbil^B ~ Mq) 


(2.31) 

(2.32) 


where we have denoted 

Xabc... = xim) ■ 

i=a,b,c,... 

In the following we neglect in comparison with lepton densities, ~ 0, which fixes Pb (if 
Hb were kept non-zero, a Legendre transform should be carried out to an ensemble with fixed 
Ub)- Charge neutrality is imposed by requiring dp/dpiq = 0. From the resulting expression, 
we can obtain lepton densities as functions of the chemical potentials as na = dp/dp^^. 

The solution to the charge neutrality condition dp/dpq = 0 reads 




Xudscb 




XeprXudscb E EcXucXdsb i=e,/i,r 


mi)Pi 


The total lepton density of flavour a is 

na = t (0) Pa + "^Xinia) K - fJ'ol ■ 
The first term accounts for the neutrino density 


(2.33) 


(2.34) 


nu^ = x(.0)Pa, (2.35) 

whereas the latter term represents the density of charged leptons of flavour a. 

It is important to note that because the three lepton asymmetries are independent of each 
other, charge conservation can be balanced by the other lepton flavours. For instance, even 
if = 0, so that p^ = 0, it is still possible to have ^ Q because tau-leptons couple to pq 
and can thereby neutralize some of the charges carried by electrons and muons. 



If, however, the different lepton densities are assumed to equilibrate (which is physically 
unlikely at T > 10 MeV [15, 16]), the situation changes. We consider this case as well, given 
that it can serve as a useful test case leading to the largest possible resonance effect [7]. For 
equilibrated lepton densities we set a G {1,2,3}, so that eq. (2.33) simplifies to 




XefirXudscb 


XefirXudscb ^cXucXdsb 


Ml 


(2.36) 


and the lepton density is = 3y;(0)/r^ + 2Xe^r(ML “ Mq)- The neutrino and 

charged lepton asymmetries read = x(0) and = 2x{'ma){^ii — Hq), respectively. 
Note that if all quark masses are made large compared with temperature, or we set W ^ 0, 
then eq. (2.36) implies that fiq = and lepton asymmetry is only carried by neutrinos. 
This is because in the absence of hadronic plasma constituents, charge neutrality would be 
violated if charged leptons were chemically equilibrated and had non-vanishing asymmetries. 


3. Practical implementation of the rate equations 
3.1. Summary of the basic setup 

For simplicity, we assume that during the period under consideration only one sterile neutrino 
is active, i.e. has an interaction rate comparable with the Hubble rate.® Then only one sterile 
neutrino flavour contributes in eq. (2.24). We denote this flavour by / = 1. With JC being 
the associated on-shell four-momentum, two rates are defined by 

|/^iaPTr[|Cp,,(/C)aJ ^ IVpTr[|Cp,,(-/C)»,] ^ 

where is the Standard Model spectral function from eq. (2.11). Then the basic equations 
to be solved, (2.21) and (2.24), take the forms {fj^ = 

fk = ^ X]{K(^l+Ma)-/fc]-Ra(^)+ K(^l-Ma)-/fc]^a(^)} > (3-2) 

a 

< = j^{[nF{Ei + ^^a) - fk]Raik) - [nF(Fli - Ma) - fk]Rtik)] , (3.3) 

with //^’s and n„’s related through eq. (2.34). The Fermi distributions in eq. (3.2), which 
are always below unity, take care of Pauli blocking. The two terms on the right-hand sides 
of these equations can be interpreted as originating from reactions involving SM “particles” 

®In the ‘TMSM”, the two heavier sterile neutrinos are assumed to have masses in the GeV range [4, 5], so 
they are not particularly “heavy” at T ~ 1 GeV. However their interaction rate peaks at higher temperatures, 
T ^ 10 — 20 GeV, cf. fig. 3 of ref. [8], and is very small at T < 1 GeV. 
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and “antiparticles”.® In the case of equilibrated active flavours, eq. (3.3) is replaced with 


nL = '^ jj^[nF{E^+fii) - fklKik) - K(£^i - - 4 ] (^)} • (3-4) 

Our goal is to solve these equations (generalized to an expanding background) in the regime 
Ml = 7.1 keV <C £^1 ~ T ~ 200 MeV <C u ~ 246 GeV , (3.5) 

where v denotes the Higgs vacuum expectation value. 


3.2. Active neutrino properties 

Given that we are deep in the Higgs phase, the Higgs doublet (j) can to a good approximation 
be replaced by its (gauge-fixed) expectation value, cf. eq. (3.5). Then the neutrino Yukawa 
couplings only appear through neutrino Dirac masses, defined as 

^ ^ . (3.6) 

The spectral function is proportional to the spectral function of an active neutrino of 
flavour a. This assignment refers to the weak interaction eigenbasis. The (retarded) propa¬ 
gator of an active neutrino is of the form [25] 

S-\-lC)=lf,a+ii {b + c + i^ , (3.7) 

and the spectral function is given by p{K,) = ImS'(/C -|- iuO'^), where u = (1, 0) is the four- 
velocity of the heat bath. The function a can to a good approximation be replaced by its 
tree-level value a = 1 [26]. The functions b and c are often called (thermal and asymmetry 
induced) matter potentials, and T can be called (up to trivial factors) a thermal width, a 
damping rate, an interaction rate, or an opacity. The function b is real, even in and odd 
in fC. The function c is real, odd in and even in JC. The imaginary part T can to a good 
approximation be evaluated at p^^ = 0 and is then even in fC. Inserting the form of eq. (3.7) 
into eq. (3.1), accounting properly for various sign conventions, and dropping terms which 
are subleading in the regime of eq. (3.5), we obtain 


R-{k) 




\M? 2£. (b + c) + (b 


R+(k) = R-(k)\ 


In o\ 


®More precisely, the rate describes a transition lepton -o- Yi and the rate Rt a transition antilepton 
Yj. Time can run in either direction. Eq. (3.2) states that sterile neutrinos can be produced from leptons 
and antileptons. Eq. (3.3) states that asymmetries between lepton and antilepton densities decrease by the 
difference of the rates felt by leptons and antileptons. 
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Let us now discuss the explicit forms of the functions appearing in eq. (3.7). In the regime of 
eq. (3.5) the imaginary part L originates at 2-loop level, and has been computed with account 
of all SM reactions in ref. [27] (there it was represented by the combination Iq ~ LiiL). It is 
of the form 

T = GlT^EiiQ, (3.9) 

where = g'^/(4\/2m,^) is the Fermi constant, and the dimensionless function Ig ~ 1 
has been tabulated for various momenta, temperatures, and lepton flavours on the web page 
related to ref. [27]. 

The function b originates at 1-loop level and was determined in ref. [26] for <C It 
can be expressed as 


b 

(p{m) 


IGGjE^ 

TTOi^ 

f d^p Up 


(27r)3 2E 


cos^Oy^ (p{0) -I- 2(/)(mQ 

+ </>( 0 ) 


77r2r4 

360 


(3.10) 

(3.11) 


where = g‘^/{A7r) and 6*^ is the weak mixing angle. In analogy with eq. (3.9) we can write 


b = GlT‘^Ei b . 


(3.12) 


Because of the l/a^, factor, the dimensionless function b is much larger than Iq, b ~ 80. It 
is plotted in fig. 1 of ref. [27] for different flavours and temperatures. 

The last ingredient is the function c, which incorporates effects from charge asymmetries. 
This function was also determined in ref. [26]. Assuming chiral equilibrium and including all 
light SM particles, we obtain 


c = 


V2G^ 2n^^+^n^^ + (i + 2 sin^^^^n^^ - - 2 sin^^^) ^ 

b^a b^a 


-sm"' 

\2 3 


E 

i=u^c 


n.-l---sm 


i=d,s,b 


(3.13) 


where the second line encodes the contribution of quarks, coming from tadpole diagrams 
mediated by the Z boson. Because the latter couples differently to up- and down-type quarks, 
the hadronic contribution contains a part that is not proportional to ^ '^i=u ds cb'^i 

and thus survives even if, as we assume, the baryon density vanishes. The quark densities 
read 


2A(.(^^ -|- 2g,Q)Xy 


nd = 


2A"c(Mb ^^Q)Xd 


(3.14) 


and correspondingly for other up- and down-type quarks, respectively. Eq. (2.31) can be used 


11 



for expressing in terms of /Xq, yielding (for —)■ 0) 

c = \/2Gp + +2sin20^)ne^- -2sin20^)^ng^ 

b^a b^a 

+2(1 - 2sin^g^) ^^^^'^^'^"^ Hq . (3.15) 

^udscb 

We again consider two cases, corresponding to those in sec. 2.4. For unequilibrated lepton 
asymmetries, the neutrino asymmetries are different, as given by eq. (2.35). Then /Xg can 
be read off from eq. (2.33). The charged lepton densities originate from the second term in 
eq. (2.34). In contrast, for equilibrated lepton asymmetries, all neutrino densities are equal, 
'iT'u — x(0) whereas Hq can be read off from eq. (2.36). 


3.3. Expanding background 


In an expanding background, the left-hand sides of eqs. (3.2), (3.3) become 

fk ~ Hkd^.)ff, , fi^ —)• + 3H)na , 


(3.16) 


where H is the Hubble rate, H = y 87re/(3mp[), and e denotes the energy density. The 
inhomogeneous term can be eliminated from the equation of motion for by integrating 
along a trajectory of redshifting momentum, 

r s(T) 1 

kr^K , (3.17) 

where s is the total entropy density, and from that for by normalizing by s, 


(3.17) 


Y (T) = 

a[ )- ■ 


(3.18) 


It is also convenient to integrate in terms of the temperature T rather than the time t. 
Denoting the final moment of integration by = 1 MeV, we get 


[nAEi + Ha) - fk^] (kr) + [M^i - Ha) “ fkj (^r) 
dln(Tjr) “ ^ GH{T)cUT) 


(3.19) 


dln(Tjr) ^ 6H{T)c1{T) > V • X 

dy„(r) _ r [MEi + Ha) - fkj Ra (^t) - K(^l - Ha) - fk^] Rt (^t) 
dln(Tjr) - ik, 3s{T)H{T)cj{T) ’ ^ ^ 

where is the speed of sound squared. Numerical values for the thermodynamic functions 
appearing in these equations (e,s,Cg) have been tabulated in ref. [28]. Note that the right- 
hand side of eq. (3.19) is even in charge conjugation, whereas that of eq. (3.20) is odd. For 
the case of equilibrated active flavours, eq. (3.20) gets replaced with (Y^ = Ya) 


dYLiT) 

dln(TjT) 


\n^{Ei + Hl) ~ fk^. Ra {krp) — — Hl) — (^ t ) 

2,s{T)H{T)cl{T) 


. (3.21) 
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3.4. Resonance contributions 


For small F the rates in eq. (3.8) resemble Dirac delta-functions, if the first factor in the 
denominator has a zero. For positive c, this can be the case with the term . Denoting 

R{T) = Ml + 2E^{h - |c|) + (6 - \c\f , (3.22) 


we can then approximate 

“ El E^ > Q 1 y > 

If c < 0, a similar term exists in R~] it can only exist in one of the terms at a time. 

In general, there are two resonances to be considered. The simplest way to see this is to 
fix T and consider T" as a function of Ei. Indeed the energy dependence of the variables 
appearing in eq. (3.22) is simple: b is linear in E^ whereas c is constant (cf. eqs. (3.10), 
(3.15)). If we write b = bE^, where 6 <C 1, then 


E=b{2 + b)El - 2Ei\c\ (1 + 6) + Ml + , 

and zeros exist if > b{2 + b)Ml. The zeros are located at 




c|(I + 6)±-^c2-6(2 + 6)M2 

6 (2 + 6 ) 


(3.24) 


(3.25) 


and the “Jacobian” reads \d^^E\ = 2 ^— 6 (2 + h)Ml at E^ = E^^. 

In practice, of course, Iq is not infinitesimally small, and the resonance is not arbitrarily 
narrow. In this situation resonance effects interfere with non-resonant contributions. One 
way to account for this in a practical numerical solution is sketched in appendix A. 


3.5. Relic density 


Once the final spectrum has been obtained through the integration of eqs. (3.19)-(3.21), 
we need to relate it to the present-day dark matter energy density. In practice we choose the 
lowest temperature of the integration to be = 1 MeV, by which time all the source terms 
have switched off, whereas Tq denotes the present-day temperature of the cosmic microwave 
background. Today, the sterile neutrinos are non-relativistic, so that the energy density 
carried by them reads 



(3.26) 


The dark matter energy density can be written as = P Hm /Pcr) 


Pdn 


= ^dmh X 


h^HTo 


X s(Tq 


(3.27) 


13 



sm^(2e) = 7 * 10‘", case a, Q, / = 1 

' 1 DM 



sin^(20) = 7 * 10‘", case e, H, / = 1 

' 1 DM 



Figure 1: Resonance locations in the different flavour channels for case (a) (left) and case (e) (right) 
[only a = e has a physical effect in these cases because only 7 ^ 0]. In each case we have set 
sin^(20) = 7 X 10“^^ and tuned the initial asymmetry to the value producing the correct dark matter 
abundance, cf. table 1. We have considered comoving momenta kj^ < 12 .at T), = 1 MeV; the 
continuous line indicates the upper edge of this range. 


where is the critical energy density and s(Tq) = 2 891/cm^ is the current entropy density. 
Making use of the known value of p^j. yields /9 j.j./[/i^s(Tq)] = 3.65 eV. Recalling that = 

0.12 according to Planck data [29], and dividing eq. (3.26) by eq. (3.27), we get 


2 Mi 

~ 0.12 X 3.65 eV 


d^kp /fcp 
( 277)3 


= 6950 X 


Ml 

7.1 keV 


X 


d3k,> 

(27r)3 n ’ 


(3.28) 


where we made use of the facts that f d3kj, ff^^/s{T) is temperature-independent at T <T^ 
and that s{T^) « 4.67T3. So, given the known , eq. (3.28) allows us to determine 


4. Numerical results 

We have integrated eqs. (3.19)“(3.21) numerically for a number of parameter values, starting 
at T = = 4 GeV where we assume = 0, and stopping at T = = 1 MeV where all 

source terms have switched off. Subsequently we determine the observables defined in sec. 3.5 
for Ml = 7.1 keV. To illustrate our results, let us focus on the following cases: 

(a) at T = only /ii^ ^ 0; equilibrated active flavours. 

(b) at T = Tmax! ^le 7^ d; non-equilibrated active flavours. 
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T/GeV 


sin\20) = 7 * 10 ", case a, Cl, / Cl,^., = 1 

^ ’ ’1 DM 


sin^(20) = 7 * 10 ", case e. Cl, / Cl„„, = 1 

^ ^ ’ ’1 DM 




Figure 2: The evolution of the lepton asymmetries for case (a) (left) and case (e) (right). The 
parameters are like in fig. 1. In case (a) Y^,Y^ grow initially, even though the source terms are 
not active yet, because charged r-leptons cannot carry their share of the asymmetry when T <C m.,. 
(y^ = Y^ is constant). In case (e) such a re-distribution is not possible and Y^ and Y^ are exactly 
conserved. The values of the initial neutrino asymmetries are given in table 1; the values of the 
corresponding lepton asymmetries follow from eqs. (2.33)-(2.35). Lepton asymmetries 

would be expected to equilibrate below T = 10 MeV [15, 16], in the region shown by a grey band, 
however the rates have switched off by then so this has no effect on sterile neutrino distributions. 


(c) at T = only ^ 0; equilibrated active flavours. 

(d) at T = only ^ 0; non-equilibrated active flavours. 

(e) only / 0 at T = only ^ 0; non-equilibrated active flavours. 

(f) only / 0 at T = only ^ 0; non-equilibrated active flavours. 

(g) only / 0 at T = only ^ 0; non-equilibrated active flavours. 

(h) only 0 at T = only hi^ 0; non-equilibrated active flavours. 

(i) only yf^ 0 at T = only ^ 0; non-equilibrated active flavours. 

(j) only yf^ 0 at T = only hi^ yf^ 0; non-equilibrated active flavours. 

Let us reiterate that in the case of equilibrated active flavours, one would have to assume 
active neutrino oscillations to proceed much faster than the processes considered in the present 
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sin\20) = 7 * 10 ", case a, = 1 sin^(26) = 7 * 10 ", case e, = 1 



Figure 3: Examples for the evolution of the right-handed neutrino distribution for case (a) (left) 
and case (e) (right), assuming that fk^{T = 4 GeV) = 0. The final temperature is = 1 MeV, and 
= krp denotes momenta at this temperature. The parameters are like in figs. 1, 2. For smallish 
most of the production takes place at the lower resonance temperature (cf. fig. 1). 

paper, which is unlikely to happen at T > 10 MeV [15, 16]. Nevertheless we display the results 
in order to allow for a comparison with ref. [7], to be performed in sec. 5. 

The initial state is parametrized by the neutrino asymmetry normalized to the entropy 
density, n^^js. The mixing angles are parametrized through 

401, ( 4 . 1 ) 

a=e,fi,T 1 

which is the combination that appears in the (inclusive) decay rate of sterile neutrinos to 
an active neutrino and a photon. We consider the value sin^(20) ~ 7 x 10“^^ mentioned in 
ref. [10] and the limits of sin^(20) ~ (2 — 20) x 10“^^ from ref. [11]. Confining effects are 
modelled through the phenomenological replacement -^ceff suggested in ref. [27]. 

(In ref. [27] it was checked that this recipe is consistent with Chiral Perturbation Theory at 
low T; unfortunately Chiral Perturbation Theory is not applicable at T > 100 MeV.) 

In fig. 1, the two resonance locations (in each channel) are shown for the cases (a) and (e). 
In fig. 2, the evolution of the densities is shown, and in fig. 3 the same is done for the 
distribution function fj^^. The ratio from eq. (3.28) is illustrated in fig. 4, whereas 

the differential shape of at T = = 1 MeV can be inferred from figs. 5 and 6 . The 

initial neutrino densities yielding the correct dark matter abundances in all cases (a)-(j) are 
summarized in table 1. It is remarkable that despite quite different asymmetries (cf. table 1), 
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Figure 4: The total energy density carried by sterile neutrinos today, normalized to the dark matter 
density, as a function of the initial lepton asymmetry, for case (a) (left) and case (e) (right). The 
couplings span the range indicated by refs. [10, 11]. 


cases (a), (b), (e), (f), (h) and (i) produce very similar spectra (cf. figs. 5, 6 ). 


5. Conclusions 

In view of the exciting (if unconsolidated) prospect of accounting for dark matter through 
7.1 keV sterile neutrinos [10, 11], the purpose of this paper has been to promote a previously 
proposed quantum field theoretic framework [7] from a qualitative towards a more quanti¬ 
tative level. In order to reach this goal, two types of “back reactions” (i.e. non-linearities) 
entering the basic equations have been derived from stated assumptions (cf. secs. 2.2, 2.3). 
The relation of the lepton densities and lepton chemical potentials entering these equations 
has been systematically worked out to leading order in small chemical potentials (cf. sec. 2.4). 
The equations have been written in a form which separately tracks three different flavours 
of non-equilibrated lepton densities (cf. eqs. (3.19), (3.20)). Finally, the equations have been 
numerically solved “as is”, without imposing further model assumptions at this stage. 

In a previous study [7], which relied otherwise on similar approximations as the present 
one, it was assumed that all active flavours are in chemical equilibrium, and that both the 
charged and the neutral leptons carry the same asymmetry, so that the total initial lepton 
asymmetry is effectively . This leads to a large coefficient c and correspondingly to 

a maximally efficient resonant contribution. In reality, as we have discussed, charged leptons 
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Figure 5: The sterile neutrino distribution function at T < = 1 MeV, normalized to the Fermi 

distribution nF(fc) = l/[exp(fc/T) + 1], for the initial lepton asymmetry producing precisely the 
observed dark matter energy density, for case (a) (left) and case (e) (right). Given that / <C up, 
sterile neutrinos are far below equilibrium despite their efficient resonant production. 

cannot carry such large asymmetries because of electric charge neutrality. Concretely, this 
implies that we need a larger active-sterile mixing angle or initial asymmetry for a comparable 
effect. For instance, for case (a), where we find the initial asymmetry = 12.25 x 10“® 

for sin^(20) = 7 x 10“^^, the analysis of ref. [7] would have produced the correct dark matter 
abundance already with ny^js = 9.05 x 10“® for sin^(20) = 7x 10“^^, or already for sin^(20) = 
1.5 X 10“^^ with s = 12.25 x 10“®. In other words, the difference between ref. [7] and 
the present work is of order unity. On the logarithmic scale of figs. 5, 6 the distributions of 
ref. [7] do however bear some similarity with ours, if considered at the same value of sin^(20). 

Our numerical results have been presented in sec. 4. The final spectra for all the cases 
considered can also be downloaded from http://www.laine.itp.unibe.ch/dmpheno/. It 
remains a theoretical challenge to confirm whether some of the pre-existing neutrino asym¬ 
metries in table 1 can indeed be produced by mechanisms such as the one described in ref. [8]. 

Despite the improvements of the present paper, it should be acknowledged that the solu¬ 
tion still contains theoretical uncertainties. The reason is that most of the sterile neutrino 
production takes place at temperatures of a few hundred MeV (cf. fig. 3), where hadronic 
effects play a significant role. In our work, hadronic effects have been handled through a phe¬ 
nomenological recipe introduced in ref. [27], which does correctly incorporate the fact that 
QCD displays a rapid but smooth crossover rather than an actual phase transition. Then 
hadronic uncertainties remain on a level of some tens of percent as discussed previously [7]. 
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Initial neutrino number density in 

units of 10 ® 

case 

sin2(2(9) = 2 X lO'^^ 

sin2(26») = 7 X lO'^^ 

sin2(26») = 20 X lO'^^ 

a 

14.14 

12.25 

10.81 

b 

19.30 

17.42 

15.65 

c 

13.47 

11.60 

10.69 

d 

19.11 

17.80 

17.15 

e 

33.45 

30.16 

27.08 

f 

102.77 

96.49 

88.72 

g 

100.56 

96.85 

94.99 

h 

82.51 

72.14 

63.60 

i 

82.34 

72.13 

63.75 

j 

30.91 

28.02 

26.65 


Table 1: Initial neutrino densities at Tmax = 4 GeV yielding the correct dark matter abundance, 
expressed as where x denotes the flavour relevant for the cases (a)-(j) (cf. sec. 4). For reasons 

of numerical reproducibility more digits have been shown than is the expected theoretical accuracy of 
our analysis (errors are expected on the 10 — 20% level, mainly from hadronic uncertainties). 

Eventually, if the sterile neutrino dark matter scenario establishes itself, many of these uncer¬ 
tainties can be reduced through lattice Monte Carlo measurements. As has been outlined in 
ref. [17] and in the present paper, lattice input is needed for the equation of state, for quark 
number susceptibilities, and for mesonic correlation functions in various quantum number 
channels. The first two ingredients would already be available but it is not clear whether 
including lattice input in some places and not in others would consistently improve on our 
results.'^ Nevertheless their gradual inclusion seems to present an interesting challenge for 
future work. 
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Appendix A. Practical treatment of resonance contribntions 

In this appendix we sketch one example for how the resonance contributions, outlined in 
sec. 3.4, can be handled in practice. 

was verified in ref. [27] that by the time Chiral Perturbation Theory is applicable, which would permit 
for an analytic treatment of hadronic effects, hadronic effects are below the 1% level. 
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In general, the resonances cannot be considered as infinitely narrow, and their treatment 
cannot be separated from that of non-resonant contributions. Consider the vicinity of a 
structure like in eq. (3.23). Numerical integrations take place on a discrete grid. Suppose 
that between two grid points, x^_i and Xj, where x can be chosen as ln(T^/T) in eq. (3.19) 

and as Ei on the right-hand side of eq. (3.20) = ^Jk‘^ + M^), we find that E has changed 

its sign; let Xq denote the location of the zero. In the case of the integral on the right-hand 
side of eq. (3.20), the correct contribution would be 


f 


dx (/>(xq 


-^q(®o) 


4>{x, 


0 ) 


arctan 


[J^'(xo)(x-Xo)]2 +/2(a;o) 
J"\Xq)\{Xq - Xi_i) 




whereas naively we could have estimated this through 
_ Xj - x^_i [ <i){Xi_i)lQ{x^_i) 


-|- arctan 


\^'{xq)\{x^ - Xq) 


(P(.x^)Iq{x^) 

\E^ix,_,) + lUx^_,) + E^ix^) + lUx,) 


(A.l) 


(A.2) 


In order to correct for the error, we may subtract 5 and add A to the result. If Xq is close to 
Xj or Xj_^, neighbouring cells need to be corrected as well. 

In the case of eq. (3.19), the integrated result has the structure 


f{Xi) = f{xi_i)+Y^ 
- f{xi_i) + Y, 2 


dx' 


- f{x')Xa{x') 


(A.3) 


a '^i—i 


^a'iXi-l) - fiXi_i)Xa'iXi-l) + “ fiXi)Xa'iXi) 


+ 


- f{xQ)Xr{Xo) 


arctan 


^qi^o) 


-|- arctan 


\^'{xp)\{x^-Xp) 

^qixp) 


where a' enumerates non-resonant terms; r is the resonant contribution; (/> = cfylg/^E'^ + Iq)', 
and X = xlqli^'^ ^q)- Assuming the non-resonant contribution to be subleading, which 
can be arranged by choosing Xj — x^_^ sufficiently small, the unknown value f{xp) can be 
estimated from 

/(Xj) - f{xp) 


arctan 


^q(*o) 


f{Xp) - /(Xj_i) 

This implies that we can write 


arctan 




/(xo)< arctan 


= /(Xj) arctan 


^q{xp) 
l-^'(a;o)|(xo -Xj_^) 

Iqixp) 


-|- arctan 




^qi^o) 


+ fixi_i) arctan 


l-^'(a;o)|(a;» -a^o) 
Iqixp) 


(A.4) 


(A.5) 
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Inserting this into eq. (A.3) we can solve for f{x^) in terms of the known including 

now the non-resonant contributions as well. 
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Figure 6 : Like fig. 5 but for the other cases [for ease of comparison, case (e) is reproduced here]. 
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